clear all;
clc;
%%
% Hochrechnung auf eine Sekunde Bestrahlung (Soll: 10^16 Photonen)
simfaktor = 1e10;
nDiv = 400;
nDivX = nDiv;
nDivY = nDiv;
nDivZ = nDiv;

%%
if(exist('filename', 'var') == 0)
    [filename, pathname] = uigetfile('.log');
    file = fopen(fullfile(pathname, filename), 'r');
    
    % Datedatei mit Kopf-Daten
    % Name (xCoord*nDivY*nDivZ+yCoord*nDivZ+zCoord) eDeposit doseDeposit
    input = fscanf(file,'%u\t%e\t%e\n',[3 inf]);
    
    % Bitte daran denken das File zu schliessen, sonst gibt es Probleme
    fclose(file);
end

%matrix3D = zeros(nDivX, nDivY, nDivZ);

entries = numel(input(1,:));
indexX = zeros(1,entries);
indexYZ = zeros(1,entries);

for n=1:entries
    nx = floor(input(1,n) / (nDivY*nDivZ));
    ny = floor((input(1,n) - nx*nDivY*nDivZ) / nDivZ);
    nz = input(1,n) - nx*nDivY*nDivZ - ny*nDivZ;
    
    %matrix3D(nx+1, ny+1, nDivZ-nz) = input(3,n);
    
    indexX(n) = nx+1;
    indexYZ(n) = ny*nDivZ+nDivZ-nz;
end
data = sparse(indexX, indexYZ, input(3,:), nDivX, nDivY*nDivZ);

%matrix3D = matrix3D .* simfaktor;

clear input;
%%
% Film festlegen mit 10fps
%film = avifile('radiation_dosimetry.avi', 'quality', 100, 'fps', 10);

% Figure Properties setzen
fig = figure('Position', [0, 0, 800, 664]);
load cmap;
colormap(cmap);
set(gca, 'XLim', [1 nDivX]);
set(gca, 'YLim', [1 nDivY]);

clim = [1e-14 1e-10]*simfaktor;

%for n=nDivZ/4-5:3*nDivZ/4+5
%for n=nDivZ/2-10:nDivZ/2+10
for n=ceil(nDivZ/2):ceil(nDivZ/2)
    
    % Kopf von der Seite anschauen
    %matrix2D = squeeze(matrix3D(n,:,:));
    matrix2D = reshape(full(data(n,:)), nDivY, nDivZ)'.*simfaktor;
    
    % Logarithmischer Plot (Matlab ausgetrickst)
    linear_axes = subplot(1,1,1);
    lin_plot = pcolor(linear_axes, matrix2D);
    
    set(gca, 'CLim', clim);
    colorbar('Yscale', 'log');
    
    set(linear_axes, 'Visible', 'off');
    set(lin_plot, 'Visible', 'off');
    
    log_axes = axes('Position', get(linear_axes, 'Position'));
    log_plot = pcolor(log_axes, log10(matrix2D));
    set(log_plot, 'EdgeColor','none');
    set(gca, 'CLim', log10(clim));
    
    axis square tight;
    
    title('Radiation Dosimetry of Patient (Gy)');
    xlabel(n)
    
    %frame = getframe(fig);
    %film = addframe(film, frame);
    print('-dpng', sprintf('output/deposition_%03d',n));
end

%film = close(film);

